home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / PMAT12 / PMATOP.PAS < prev    next >
Pascal/Delphi Source File  |  1993-01-17  |  48KB  |  1,555 lines
  1. {
  2. P-Mat - A Turbo Pascal program for Recursive Matrix Algebra
  3.  
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/30/93
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9.  
  10.  
  11. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  12. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  13. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  14. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  15. FROM USE OF THIS PROGRAM.
  16. }
  17.  
  18. Unit pmatop;
  19.  
  20. Interface
  21. Uses pmat;
  22.  
  23. Function MSort( a: vmatrixptr; col, order: integer ): vmatrixptr;
  24. Function Mexp( ROp: vmatrixptr ): vmatrixptr;
  25. Function Mabs( ROp: vmatrixptr ): vmatrixptr;
  26. Function Mlog( ROp: vmatrixptr ): vmatrixptr;
  27. Function Msqrt( ROp: vmatrixptr ): vmatrixptr;
  28.  
  29. Function Trace( ROp: vmatrixptr ): double;
  30. Function Helm( n: integer ): vmatrixptr;
  31. Function Index( lower, upper, step: integer ): vmatrixptr;
  32. Function Kron( a, b: vmatrixptr ): vmatrixptr;
  33. Function Det( Datain: vmatrixptr ): double;
  34.  
  35. Function Cholesky( ROp: vmatrixptr ): vmatrixptr;
  36. Procedure QR( Var ROp, Q, R: vmatrixptr; makeq: boolean );
  37. Function Svd( Var A, U, S, V: vmatrixptr; 
  38.              makeu, makev: boolean ) : integer;
  39. Function Ginv( ROp: vmatrixptr ): vmatrixptr;
  40. Function Fft( ROp: vmatrixptr; INZEE: boolean ): vmatrixptr;
  41.  
  42. Function Vec( ROp: vmatrixptr ): vmatrixptr;
  43. Function Diag( ROp: vmatrixptr ): vmatrixptr;
  44. Function Shape( ROp: vmatrixptr; nrows: integer ): vmatrixptr;
  45.  
  46. Type  margins = (ALL,ROWS,COLUMNS);
  47. Function Sum( ROp: vmatrixptr; marg: margins ) : vmatrixptr;
  48. Function Sumsq( ROp: vmatrixptr; marg: margins ): vmatrixptr;
  49. Function Cusum( ROp: vmatrixptr ): vmatrixptr;
  50. Function Mmax( ROp: vmatrixptr; marg: margins ): vmatrixptr;
  51. Function Mmin( ROp: vmatrixptr; marg: margins ): vmatrixptr;
  52.  
  53. Procedure Crow( Var ROp: vmatrixptr; row: integer; c: double );
  54. Procedure Srow( Var ROp: vmatrixptr; row1, row2: integer );
  55. Procedure Lrow( Var ROp: vmatrixptr; row1, row2: integer; c: double );
  56. Procedure Ccol( Var ROp: vmatrixptr; col: integer; c: double );
  57. Procedure Scol( Var ROp: vmatrixptr; col1, col2: integer );
  58. Procedure Lcol( Var ROp: vmatrixptr; col1, col2: integer; c: double );
  59.  
  60.  
  61. (************************************************************)
  62. Implementation
  63.  
  64. Procedure heapsort( Var a: vmatrixptr );
  65. Var
  66.   n, s0, s1, s2, s3, s4, s5, s4m1, s5m1 : integer;
  67.   t1, ts : double;
  68. Begin
  69.     { Shell-Metzger sort, PC-World, May 1986 }
  70.     n := a^.r;
  71.     s0 := n;
  72.     s1 := n;
  73.     s1 := s1 Div 2;
  74.     While s1 > 0 Do Begin
  75.         s2 := s0 - s1;
  76.         For s3 := 1 To s2 Do Begin 
  77.             s4 := s3;
  78.             While s4 > 0 Do Begin 
  79.                 s5 := s4 + s1;
  80.                 s4m1 := s4;
  81.                 s5m1 := s5;
  82.                 If a^.m( s4m1, 1 ) > a^.m( s5m1, 1 ) Then Begin 
  83.                     t1 := a^.m( s4m1, 1 );
  84.                     a^.mm( s4m1, 1 )^ := a^.m( s5m1, 1 );
  85.                     a^.mm( s5m1, 1 )^ := t1;
  86.                     ts := a^.m( s4m1, 2 );
  87.                     a^.mm( s4m1, 2 )^ := a^.m( s5m1, 2 );
  88.                     a^.mm( s5m1, 2 )^ := ts;
  89.                     s4 := s4 - s1;
  90.                 End
  91.                 Else Begin 
  92.                     s4 := 0;
  93.                 End;
  94.             End;
  95.         End;
  96.         s1 := s1 Div 2;
  97.     End;
  98. End; { end heapsort }
  99.  
  100.  
  101. Function MSort{(a :vmatrixptr; col, order: integer): vmatrixptr;};
  102. Var
  103.   sorted, temp : vmatrixptr;
  104.   i,j,t : integer;
  105. Begin
  106.     a^.Garbage;
  107.     new( sorted, makematrix( a^.r, a^.c ) );
  108.     If order <> 0 Then Begin 
  109.         { sort column col, then reorder other rows }
  110.         new( temp, makematrix( a^.r, 2 ) );
  111.         For i := 1 To a^.r Do Begin 
  112.             temp^.mm( i, 1 )^ := a^.m( i, col );
  113.             temp^.mm( i, 2 )^ := i ;
  114.         End;
  115.         heapsort( temp );
  116.         For i := 1 To a^.r Do Begin 
  117.             For j := 1 To a^.c Do Begin 
  118.                 t := Trunc( temp^.m( i, 2 ) );
  119.                 sorted^.mm( i, j )^ := a^.m( t, j );
  120.             End;
  121.         End;
  122.     End
  123.     Else Begin 
  124.         { sort row col, then reorder other columns }
  125.         new( temp, makematrix( a^.c, 2 ) );
  126.         For i := 1 To a^.c Do Begin 
  127.             temp^.mm( i, 1 )^ := a^.m( col, i );
  128.             temp^.mm( i, 2 )^ := i ;
  129.         End;
  130.         heapsort( temp );
  131.         For i := 1 To a^.c Do Begin 
  132.             For j := 1 To a^.r Do Begin  
  133.                 t := Trunc( temp^.m( i, 2 ) );
  134.                 sorted^.mm( j, i )^ := a^.m( j, t );
  135.             End;
  136.         End;
  137.     End;
  138.     Dispatch^.Push( sorted );
  139.     dispose( temp, killVmatrix );
  140.     MSort := Dispatch^.ReturnMat;
  141. End;
  142. Function Mexp{(ROp :vmatrixptr): vmatrixptr;};
  143. Var
  144.   temp : vmatrixptr;
  145.   i,j  : integer;
  146. Begin
  147.     {  take exponent of all elements  }
  148.     ROp^.Garbage;
  149.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  150.     For i := 1 To  ROp^.r Do Begin 
  151.         For j := 1 To  ROp^.c Do
  152.             temp^.mm( i, j )^ := exp( ROp^.m( i, j ) );
  153.     End;
  154.     Dispatch^.Push( temp );
  155.     Mexp := Dispatch^.ReturnMat;
  156. End;
  157.  
  158. Function Mabs{(ROp: vmatrixptr): vmatrixptr;};
  159. Var
  160.   temp : vmatrixptr;
  161.   i,j  : integer;
  162. Begin
  163.     {  take absolute values of all elements  }
  164.     ROp^.Garbage;
  165.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  166.     For i := 1 To  ROp^.r Do Begin    
  167.         For j := 1 To  ROp^.c Do
  168.             temp^.mm( i, j )^ := abs( ROp^.m( i, j ) );
  169.     End;
  170.     Dispatch^.Push( temp );
  171.     Mabs := Dispatch^.ReturnMat;
  172. End;
  173.  
  174. Function Mlog{(ROp :vmatrixptr): vmatrixptr;};
  175. Var
  176.   temp : vmatrixptr;
  177.   i,j  : integer;
  178.   R    : double;
  179. Begin
  180.     {  take logarithm of all elements  }
  181.     ROp^.Garbage;
  182.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  183.     For i := 1 To  ROp^.r Do Begin    
  184.         For j := 1 To  ROp^.c Do Begin 
  185.             R := ROp^.m( i, j );
  186.             If R <= 0.0 Then
  187.                 Dispatch^.Nrerror( 'Mlog: zero or negative argument to log' )
  188.             Else temp^.mm( i, j )^ := ln( R );
  189.         End;
  190.     End;
  191.     Dispatch^.Push( temp );
  192.     Mlog := Dispatch^.ReturnMat;
  193. End;
  194.  
  195. Function Msqrt{(ROp: vmatrixptr): vmatrixptr;};
  196. Var
  197.   temp : vmatrixptr;
  198.   i,j  : integer;
  199.   R    : double;
  200. Begin
  201.     {  take sqrt of all elements  }
  202.     ROp^.Garbage;
  203.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  204.     For i := 1 To  ROp^.r Do Begin 
  205.         For j := 1 To  ROp^.c Do Begin 
  206.             R := ROp^.m( i, j );
  207.             If R < 0 Then
  208.                 Dispatch^.Nrerror( 'Msqrt: zero or negative argument to sqrt' )
  209.             Else temp^.mm( i, j )^ := sqrt( R );
  210.         End;
  211.     End;
  212.     Dispatch^.Push( temp );
  213.     Msqrt := Dispatch^.ReturnMat;
  214. End;
  215.  
  216. Function Trace{(ROp: vmatrixptr): double;};
  217. Var
  218.   i : integer;
  219.   t : double;
  220. Begin
  221.     ROp^.Garbage;
  222.     t := 0;
  223.     If ROp^.r <> ROp^.c Then
  224.         Dispatch^.Nrerror( ' Trace: matrix not square in trace' );
  225.     For i := 1 To  ROp^.r Do t := t + ROp^.m( i, i );
  226.     trace := t;
  227. End;
  228.  
  229. Function Helm {(n: integer): vmatrixptr;};
  230. Var
  231.   i, j : integer;
  232.   d, den : double;
  233.   temp : vmatrixptr;
  234. Begin
  235.     new( temp, makematrix( n, n ) );
  236.     For i := 1  To n Do Begin    
  237.         For j := 1 To  n Do Begin    
  238.             d := 1.0 / sqrt( 0.0 + n );
  239.             den := d;
  240.             If j > 1 Then den := 1.0 / sqrt( 0.0 + j * (j - 1) );
  241.             If (j > 1) And  (j < i) Then d := 0;
  242.             If (j > 1) And  (j = i) Then d := - den * ( 0.0 + (j - 1));
  243.             If (j > 1) And  (j > i) Then d := den;
  244.             temp^.mm( i, j )^ := d;
  245.         End;
  246.     End;
  247.     Dispatch^.Push( temp );
  248.     helm := Dispatch^.ReturnMat;
  249. End;
  250.  
  251. Function Index{( lower,  upper, step : integer): vmatrixptr;};
  252. Var
  253.   cnter, i : integer;
  254.   temp : vmatrixptr;
  255. Begin
  256.     If step = 0 Then Dispatch^.Nrerror( 'Index: step is zero' );
  257.     cnter := 0;
  258.     i := lower;
  259.     If lower < upper Then Begin 
  260.         If step < 0 Then
  261.             Dispatch^.Nrerror( 'Index: trying to step from lower to upper with negative step size' );
  262.         While i <= upper Do Begin 
  263.             cnter := cnter + 1;
  264.             i := i + step;
  265.         End;
  266.     End
  267.     Else Begin 
  268.         If step > 0 Then
  269.             Dispatch^.Nrerror( 'Index: trying to step from upper to lower with positive step size' );
  270.         While i >= upper Do Begin 
  271.             cnter := cnter + 1;
  272.             i := i + step;
  273.         End;
  274.     End;
  275.     If cnter = 0 Then
  276.         Dispatch^.Nrerror( 'Index: trying to allocate a matrix with zero rows' );
  277.     
  278.     new( temp, makematrix( cnter, 1 ) );
  279.     
  280.     cnter := 1;
  281.     i := lower;
  282.     If lower < upper Then
  283.         While i <= upper Do Begin 
  284.             temp^.mm( cnter, 1 )^ := i;
  285.             i := i + step;
  286.             cnter := cnter + 1;
  287.         End
  288.     Else 
  289.         While i >= upper Do Begin 
  290.             temp^.mm( cnter, 1 )^ := i;
  291.             i := i + step;
  292.             cnter := cnter + 1;
  293.         End;
  294.     
  295.     Dispatch^.Push( temp );
  296.     Index := Dispatch^.ReturnMat;
  297. End;
  298.  
  299. Function Kron{(a,b:vmatrixptr): vmatrixptr;};
  300. Var
  301.   cc,cr,i,j,k,l : integer;
  302.   c : vmatrixptr;
  303. Begin
  304.     { Kroniker's product }
  305.     a^.Garbage;
  306.     b^.Garbage;
  307.     
  308.     cr := a^.r * b^.r;
  309.     cc := a^.c * b^.c;
  310.     new( c, makematrix( cr, cc ) );
  311.     
  312.     For i := 1 To  a^.r Do Begin 
  313.         For j := 1 To  a^.c Do Begin 
  314.             For k := 1 To  b^.r Do Begin 
  315.                 For l := 1 To  b^.c Do Begin 
  316.                     c^.mm( (i - 1) * b^.r + k, (j - 1) * b^.c + l )^ := a^.m( i, j ) * b^.m( k, l );
  317.                 End;
  318.             End;
  319.         End;
  320.     End;
  321.     Dispatch^.Push( c );
  322.     Kron := Dispatch^.ReturnMat;
  323. End;     (* kron *)
  324.  
  325. { private function in unit: used by determinant }
  326. Procedure Pivot( Var Data: vmatrixptr; RefRow: integer; 
  327.     Var Determ: double; Var DetEqualsZero: boolean );
  328. Var
  329.   NewRow, i: integer;
  330.   temp : double;
  331. Begin
  332.     DetEqualsZero := TRUE;
  333.     NewRow := RefRow;
  334.     While DetEqualsZero And  (NewRow < Data^.r) Do Begin 
  335.         NewRow := NewRow + 1;
  336.         If abs( Data^.m( NewRow, RefRow ) ) > 1.0E - 8 Then Begin 
  337.             For i := 1 To  Data^.r Do Begin 
  338.                 temp := Data^.m( NewRow, i );
  339.                 Data^.mm( NewRow, i )^ := Data^.m( RefRow, i );
  340.                 Data^.mm( RefRow, i )^ := temp;
  341.             End;
  342.             DetEqualsZero := FALSE;
  343.             Determ := - Determ;        (* row swap adjustment to det *)
  344.         End;
  345.     End;
  346. End;
  347.  
  348. Function Det{(Datain : vmatrixptr): double;};
  349. Var
  350.   Determ : double;
  351.   data   : vmatrixptr;
  352.   coef   : extended;
  353.   Row, RefRow, Dimen, i : integer;
  354.   DetEqualsZero : boolean;
  355. Begin
  356.     Datain^.Garbage;
  357.     Determ := 1;
  358.     If Datain^.r = Datain^.c Then Begin 
  359.         Dispatch^.Inclevel;
  360.         Dimen := Datain^.r;
  361.         new( Data, makematrix( Dimen, Dimen ) );
  362.         
  363.         Data := matequals( Data, Datain );
  364.         Row := 0;
  365.         RefRow := 0;
  366.         DetEqualsZero := FALSE;
  367.         
  368.         While (Not DetEqualsZero) And (RefRow < Dimen - 1) Do Begin 
  369.             RefRow := RefRow + 1;
  370.             If abs( Data^.m( RefRow, RefRow ) ) < 1.0E - 8 Then
  371.                 Pivot( Data, RefRow, Determ, DetEqualsZero );
  372.             If (Not DetEqualsZero) Then
  373.                 For Row := RefRow + 1 To  Dimen Do
  374.                     If abs( Data^.m( Row, RefRow ) ) > 1.0E - 8 Then Begin 
  375.                         Coef := - Data^.m( Row, RefRow ) / Data^.m( RefRow, RefRow );
  376.                         For i := 1 To  Dimen Do
  377.                             Data^.mm( Row, i )^ := Data^.m( Row, i ) +
  378.                                 (Coef * Data^.m( RefRow, i ));
  379.                     End;
  380.             Determ := Determ * Data^.m( RefRow, RefRow );
  381.         End;
  382.         If DetEqualsZero Then Determ := 0
  383.         Else Determ := Determ * Data^.m( Dimen, Dimen );
  384.         dispose( Data, killvmatrix );
  385.         Dispatch^.Declevel;
  386.     End
  387.     Else Dispatch^.Nrerror( 'Det: Matrix is not square\n' );
  388.     det := Determ;
  389. End;
  390.  
  391. {  / --------------- cholesky decomposition ----------------------  }
  392. {  /  ROp is symmetric, returns upper triangular S where ROp := S'S  }
  393. {  / -------------------------------------------------------------  }
  394.  
  395. Function Cholesky{(ROp: vmatrixptr): vmatrixptr;};
  396. Var
  397.   nr, i, j, k  : integer;
  398.   temp : vmatrixptr;
  399.   sum : double;
  400. Begin
  401.     nr := ROp^.r;
  402.     If ROp^.r <> ROp^.c Then
  403.         Dispatch^.Nrerror( 'Cholesky: Input matrix not square' );
  404.     
  405.     ROp^.Garbage;
  406.     For i := 1 To  nr Do Begin    
  407.         For j := 1 To i - 1 Do
  408.             If abs( ROp^.m( i, j ) - ROp^.m( j, i ) ) > 1.0E - 7 Then
  409.                 Dispatch^.Nrerror( ' Cholesky: Input matrix is not symmetric' );
  410.     End;
  411.     new( temp, makematrix( nr, nr ) );
  412.     temp := matequals( temp, Fill( nr, nr, 0.0 ) );
  413.     For i := 1 To  nr Do Begin    
  414.         For j := i To  nr Do Begin    
  415.             sum := 0.0;
  416.             For k := 1 To i - 1 Do
  417.                 sum := sum + temp^.m( k, i ) * temp^.m( k, j );
  418.             If i = j Then
  419.                 If ROp^.m( i, j ) < sum Then
  420.                     Dispatch^.Nrerror( ' Cholesky: negative pivot' )
  421.                 Else
  422.                     temp^.mm( i, j )^ := sqrt( ROp^.m( i, j ) - sum )
  423.             Else
  424.                 If abs( temp^.m( i, i ) ) < 1.0E - 7 Then
  425.                     Dispatch^.Nrerror( ' Cholesky: zero or negative pivot' )
  426.             Else
  427.                 temp^.mm( i, j )^ := (ROp^.m( i, j ) - sum) / temp^.m( i, i );
  428.         End;
  429.     End;
  430.     Dispatch^.Push( temp );
  431.     Cholesky := Dispatch^.ReturnMat;
  432. End;
  433.  
  434. Procedure QR{(var ROp, Q, R: vmatrixptr; makeq: boolean);};
  435. Var
  436. rr,c,j,i,k: integer;
  437. s, sigma, beta, sum  : double;
  438. Begin
  439.     ROp^.Garbage;
  440.     Q^.Garbage;
  441.     R^.Garbage;
  442.     
  443.     Dispatch^.Inclevel;
  444.     rr := ROp^.r;
  445.     c := ROp^.c;
  446.     
  447.     Q := matequals( Q, ROp );
  448.     R := matequals( R, Fill( c, c, 0.0 ) );
  449.     
  450.     For j := 1 To  c Do Begin    
  451.         sigma := 0.0;
  452.         For i := j To  rr Do
  453.             sigma := sigma + Q^.m( i, j ) * Q^.m( i, j );
  454.         If abs( sigma ) <= 1.0e - 10 Then
  455.             Dispatch^.Nrerror( ' QR: singular X' );
  456.         If Q^.m( j, j ) < 0 Then s := - sqrt( sigma )
  457.         Else s := sqrt( sigma );
  458.         R^.mm( j, j )^ := s;
  459.         beta := 1.0 / (s * Q^.m( j, j ) - sigma);
  460.         Q^.mm( j, j )^ := Q^.m( j, j ) - s;
  461.         For k := j + 1 To  c Do Begin    
  462.             sum := 0.0;
  463.             For i := j To  rr Do
  464.                 sum := sum + Q^.m( i, j ) * Q^.m( i, k );
  465.             sum := sum * beta;
  466.             For i := j To  rr Do
  467.                 Q^.mm( i, k )^ := Q^.m( i, k ) + Q^.m( i, j ) * sum;
  468.             R^.mm( j, k )^ := Q^.m( j, k );
  469.         End;
  470.     End;
  471.     
  472.     If makeq Then Q := matequals( Q, Mult( ROp, Inv( R ) ) );
  473.     Dispatch^.Declevel;
  474. End;
  475.  
  476. {  / ------------------- Singular Valued Decomposition ----------  }
  477. {  /  from EISPACK SVD  }
  478. {  / ------------------------------------------------------------  }
  479. Function safehypot( a1, b1: extended ): extended;
  480. { function to get hypotenuse numbers a1 and b1  }
  481. { where a:=ln(a1) and b:=ln(b1)  }
  482. Var
  483.   del, sum, a, b : extended;
  484. Begin
  485.     a := 2.0 * ln( abs( a1 ) );
  486.     b := 2.0 * ln( abs( b1 ) );
  487.     If a > b Then del := a 
  488.     Else del := b;
  489.     { add a1^2 + b1^2, but in logarithms  }
  490.     sum := ln( exp( a - del ) + exp( b - del ) ) + del;
  491.     safehypot := (exp( 0.5 * sum ));
  492. End;
  493.  
  494. Procedure svdtemp( a: vmatrixptr; 
  495.                   Var w: vmatrixptr; 
  496.                   matu: boolean; Var u: vmatrixptr; 
  497.                   matv: boolean; Var v: vmatrixptr; 
  498.                   Var ierr: integer; Var rv1: vmatrixptr );
  499. Var
  500.    m, n, i, j, k, l, ii, i1, kk, k1, ll, l1, its, mn : integer;
  501.    c, f, g, h, s, x, y, z, eps, scale, machep, fss : extended;
  502. {
  503.   apologies to the purists, but I wanted to preserve the
  504.   FORTRAN style here because SVD is a good use of goto's. It
  505.   will also help you check it if you go to the original code.
  506. }
  507. Label l190, l210, l270, l290, l360, l390, l410,
  508.       l430, l460, l475, l490, l510, l520, l540,
  509.       l565, l600, l650, l700, l1000;
  510. Begin
  511.     Dispatch^.Inclevel;
  512.     a^.Garbage;
  513.     w^.Garbage;
  514.     u^.Garbage;
  515.     v^.Garbage;
  516.     rv1^.Garbage;
  517.     
  518.     m := a^.r;
  519.     n := a^.c;
  520.     
  521.     { vmatrixptr  a(' a' ,m,n),w(' w' ,n,1),  }
  522.     {             u(' u' ,m,m),v(' v' ,n,n),rv1(' rv1' ,n,1);  }
  523.     {        boolean matu,matv;  }
  524.     
  525.     {        ********** Machep is a machine dependent parameter specifying  }
  526.     {        the relative precision of floating point aritmetic.  }
  527.     {        for pc doubles, range is 1.6^-308 to 1.6^308  }
  528.     {        **************  }
  529.     
  530.     machep := exp( - 308.0 * ln( 1.6 ) );
  531.     
  532.     ierr := 0;
  533.     
  534.     u := matequals( u, a );
  535.     {   /      ********householder reduction to bi diagonal form ********  }
  536.     g := 0.0;
  537.     scale := 0.0;
  538.     x := 0.0;
  539.     
  540.     For i := 1 To  n Do Begin    
  541.         l := i + 1;
  542.         rv1^.mm( i, 1 )^ := scale * g;
  543.         g := 0.0;
  544.         s := 0.0;
  545.         scale := 0.0;
  546.         If i > m Then Goto l210;
  547.         
  548.         For k := i To m Do scale := scale + abs( u^.m( k, i ) );
  549.         
  550.         If scale = 0.0 Then Goto l210;
  551.         
  552.         For k := i To m Do Begin 
  553.             u^.mm( k, i )^ := u^.m( k, i ) / scale;
  554.             s := s + u^.m( k, i ) * u^.m( k, i );
  555.         End;
  556.         
  557.         f := u^.m( i, i );
  558.         {g := -((f >= 0.0) ? abs(sqrt(s)) : - abs(sqrt(s)));}
  559.         If f >= 0.0 Then g := - abs( sqrt( s ) )
  560.         Else g := abs( sqrt( s ) );
  561.         h := f * g - s;
  562.         u^.mm( i, i )^ := f - g;
  563.         If i = n Then Goto l190;
  564.         
  565.         For j := l To n Do Begin 
  566.             s := 0.0;
  567.             For k := i To m Do s := s + u^.m( k, i ) * u^.m( k, j );
  568.             f := s / h;
  569.             For k := i To m Do u^.mm( k, j )^ := u^.m( k, j ) + f * u^.m( k, i );
  570.         End;
  571.         l190 : 
  572.             For k := i To m Do u^.mm( k, i )^ := scale * u^.m( k, i );
  573.         
  574.         l210 : w^.mm( i, 1 )^ := scale * g;
  575.         g := 0.0;
  576.         s := 0.0;
  577.         scale := 0.0;
  578.         If (i > m) Or (i = n) Then Goto l290;
  579.         
  580.         For k := l To n Do scale := scale + abs( u^.m( i, k ) );
  581.         
  582.         If scale = 0.0 Then Goto l290;
  583.         
  584.         For k := l To n Do Begin 
  585.             u^.mm( i, k )^ := u^.m( i, k ) / scale;
  586.             s := s + u^.m( i, k ) * u^.m( i, k );
  587.         End;
  588.         f := u^.m( i, l );
  589.         {g := -((f >= 0.0) ? abs(sqrt(s)) : - abs(sqrt(s)));}
  590.         If f >= 0.0 Then g := - abs( sqrt( s ) )
  591.         Else g := abs( sqrt( s ) );
  592.         h := f * g - s;
  593.         u^.mm( i, l )^ := f - g;
  594.         
  595.         For k := l To n Do rv1^.mm( k, 1 )^ := u^.m( i, k ) / h;
  596.         
  597.         If i = m Then Goto l270;
  598.         
  599.         For j := l To m Do Begin 
  600.             s := 0.0;
  601.             For k := l To n Do s := s + u^.m( j, k ) * u^.m( i, k );
  602.             For k := l To n Do u^.mm( j, k )^ := u^.m( j, k ) + s * rv1^.m( k, 1 );
  603.         End;
  604.         
  605.         l270 : 
  606.             For k := l To n Do u^.mm( i, k )^ := scale * u^.m( i, k );
  607.         
  608.         l290 :
  609.         {x := (x > abs(w^.m(i, 1)) + abs(rv1^.m(i, 1))) ? x :
  610.             abs(w^.m(i, 1)) + abs(rv1^.m(i, 1)); }
  611.             If x <= abs( w^.m( i, 1 ) ) + abs( rv1^.m( i, 1 ) ) Then
  612.                 x := abs( w^.m( i, 1 ) ) + abs( rv1^.m( i, 1 ) );
  613.     End;
  614.     {       ********** accumulation of right-hand transformations ********  }
  615.     
  616.     If Not matv Then Goto l410;
  617.     
  618.     For ii := 1 To n Do Begin 
  619.         i := n + 1 - ii;
  620.         If i = n Then Goto l390;
  621.         If g = 0.0 Then Goto l360;
  622.         
  623.         { double division avoids possible underflow  }
  624.         For j := l To n Do v^.mm( j, i )^ := (u^.m( i, j ) / u^.m( i, l )) / g;
  625.         
  626.         For j := l To n Do Begin 
  627.             s := 0.0;
  628.             For k := l To n Do s := s + u^.m( i, k ) * v^.m( k, j );
  629.             For k := l To n Do v^.mm( k, j )^ := v^.m( k, j ) + s * v^.m( k, i );
  630.         End;
  631.         
  632.         l360 : 
  633.             For j := l To  n Do Begin 
  634.                 v^.mm( i, j )^ := 0.0;
  635.                 v^.mm( j, i )^ := 0.0;
  636.             End;
  637.         
  638.         l390 : v^.mm( i, i )^ := 1.0;
  639.         g := rv1^.m( i, 1 );
  640.         l := i;
  641.     End;
  642.     {*************accumulation of left-hand transformations  }
  643.     l410 : 
  644.         If Not matu Then Goto l510;
  645.     {   /        ************* for i:= min(m,n) step -1 until 1 do ******  }
  646.     mn := n;
  647.     If m < n Then mn := m;
  648.     
  649.     For ii := 1 To mn Do Begin 
  650.         i := mn + 1 - ii;
  651.         l := i + 1;
  652.         g := w^.m( i, 1 );
  653.         If i = n Then Goto l430;
  654.         
  655.         For j := l To n Do u^.mm( i, j )^ := 0.0;
  656.         
  657.         l430 : 
  658.             If g = 0.0 Then Goto l475;
  659.         If i = mn Then Goto l460;
  660.         
  661.         For j := l To n Do Begin 
  662.             s := 0.0;
  663.             
  664.             For k := l To m Do s := s + u^.m( k, i ) * u^.m( k, j );
  665.             { double division avoids possible underflow  }
  666.             f := (s / u^.m( i, i )) / g;
  667.             For k := i To m Do u^.mm( k, j )^ := u^.m( k, j ) + f * u^.m( k, i );
  668.         End;
  669.         
  670.         l460 : 
  671.             For j := i To m Do u^.mm( j, i )^ := u^.m( j, i ) / g;
  672.         
  673.         Goto l490;
  674.         
  675.         l475 : 
  676.             For j := i To m Do u^.mm( j, i )^ := 0.0;
  677.         
  678.         l490 : u^.mm( i, i )^ := u^.m( i, i ) + 1.0;
  679.     End;
  680.     
  681.     {       ********** diagonalization of the bidiagonal form ***********  }
  682.     
  683.     l510 : eps := machep * x;
  684.     {  ********** for k:=n step -1 until 1 do ***********************  }
  685.     For kk := 1 To n Do Begin 
  686.         k1 := n - kk;
  687.         k := k1 + 1;
  688.         its := 0;
  689.         {  ********** test for splitting .**********  }
  690.         {  for l:=k step -1 until 1 do  }
  691.         l520 : 
  692.             For ll := 1 To k Do Begin 
  693.                 l1 := k - ll;
  694.                 l := l1 + 1;
  695.                 If abs( rv1^.m( l, 1 ) ) <= eps Then Goto l565;
  696.                 {   /   rv1(1) is always zero, so there is no exit  }
  697.                 {   /   through the bottom of the loop *********  }
  698.                 If abs( w^.m( l1, 1 ) ) <= eps Then Goto l540;
  699.             End;
  700.         {  ************* cancellation of rv1(l) if l greater than 1******  }
  701.         l540 : c := 0.0;
  702.         s := 1.0;
  703.         
  704.         For i := l To k Do Begin 
  705.             f := s * rv1^.m( i, 1 );
  706.             rv1^.mm( i, 1 )^ := c * rv1^.m( i, 1 );
  707.             If abs( f ) <= eps Then Goto l565;
  708.             g := w^.m( i, 1 );
  709.             h := safehypot( f, g );
  710.             w^.mm( i, 1 )^ := h;
  711.             c := g / h;
  712.             s := - f / h;
  713.             If matu Then Begin 
  714.                 { go to 560 }
  715.                 For j := 1 To m Do Begin 
  716.                     y := u^.m( j, l1 );
  717.                     z := u^.m( j, i );
  718.                     u^.mm( j, l1 )^ := y * c + z * s;
  719.                     u^.mm( j, i )^ := - y * s + z * c;
  720.                 End;
  721.             End;
  722.         End;
  723.         {  ************** test for convergence ********************  }
  724.         l565 : z := w^.m( k, 1 );
  725.         If l = k Then Goto l650;
  726.         {  *************** if shift from bottom 2 by 2 minor *******  }
  727.         If its = 50 Then Goto l1000;
  728.         its := its + 1;
  729.         x := w^.m( l, 1 );
  730.         y := w^.m( k1, 1 );
  731.         g := rv1^.m( k1, 1 );
  732.         h := rv1^.m( k, 1 );
  733.         f := ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
  734.         g := safehypot( f, 1.0 );
  735.         {fss :=(f >= 0.0) ? abs(g) : - abs(g);}
  736.         If f >= 0.0 Then fss := abs( g )
  737.         Else fss := - abs( g );
  738.         
  739.         f := ((x - z) * (x + z) + h * (y / (f + fss) - h)) / x;
  740.         {  **************** next qr transformation ***************  }
  741.         c := 1.0;
  742.         s := 1.0;
  743.         
  744.         For i1 := l To k1 Do Begin 
  745.             i := i1 + 1;
  746.             g := rv1^.m( i, 1 );
  747.             y := w^.m( i, 1 );
  748.             h := s * g;
  749.             g := c * g;
  750.             z := safehypot( f, h );
  751.             rv1^.mm( i1, 1 )^ := z;
  752.             c := f / z;
  753.             s := h / z;
  754.             f := x * c + g * s;
  755.             g := - x * s + g * c;
  756.             h := y * s;
  757.             y := y * c;
  758.             If matv Then Begin 
  759.                 { goto l575; }
  760.                 For j := 1 To n Do Begin 
  761.                     x := v^.m( j, i1 );
  762.                     z := v^.m( j, i );
  763.                     v^.mm( j, i1 )^ := x * c + z * s;
  764.                     v^.mm( j, i )^ := - x * s + z * c;
  765.                 End;
  766.             End;
  767.             z := safehypot( f, h );
  768.             w^.mm( i1, 1 )^ := z;
  769.             { ************ rotation can be arbitrary if z is zero *******  }
  770.             If z <> 0.0 Then Begin 
  771.                 { goto l580; }
  772.                 c := f / z;
  773.                 s := h / z;
  774.             End;
  775.             { l580:  }
  776.             f := c * g + s * y;
  777.             x := - s * g + c * y;
  778.             If matu Then Begin 
  779.                 {//goto l600; }
  780.                 For j := 1 To m Do Begin 
  781.                     y := u^.m( j, i1 );
  782.                     z := u^.m( j, i );
  783.                     u^.mm( j, i1 )^ := y * c + z * s;
  784.                     u^.mm( j, i )^ := - y * s + z * c;
  785.                 End;
  786.             End;
  787.             
  788.             l600 : x := x; 
  789.         End;                           {//continue}
  790.         
  791.         rv1^.mm( l, 1 )^ := 0.0;
  792.         rv1^.mm( k, 1 )^ := f;
  793.         w^.mm( k, 1 )^ := x;
  794.         Goto l520;
  795.         {   /        ************** convergence **************  }
  796.         l650 : 
  797.             If z >= 0.0 Then Goto l700;
  798.         {   /        *************  w(k) is made non-negative *********  }
  799.         w^.mm( k, 1 )^ := - z;
  800.         If Not matv Then Goto l700;
  801.         
  802.         For j := 1 To n Do v^.mm( j, k )^ := - v^.m( j, k );
  803.         
  804.         l700 : x := x; 
  805.     End;
  806.     
  807.     
  808.     ierr := 0;
  809.     Dispatch^.Declevel;
  810.     exit;
  811.     { ***************** set error -- no convergence to a  }
  812.     { singular value after 30 iterations  }
  813.     l1000 : ierr := k;
  814.     Dispatch^.Declevel;
  815.     exit;
  816. End;
  817.  
  818.  
  819. Function Svd{(var A, U, S, V : vmatrixptr; makeu, makev: boolean;) : integer;};
  820. Var
  821.   aa, uu, ss, vv, rv1 : vmatrixptr;
  822.   ierr : integer;
  823. Begin
  824.     Dispatch^.Inclevel;
  825.     
  826.     A^.Garbage;
  827.     ierr := 0;
  828.     
  829.     new( aa, makematrix( 1, 1 ) );
  830.     new( uu, makematrix( 1, 1 ) );
  831.     new( ss, makematrix( 1, 1 ) );
  832.     new( vv, makematrix( 1, 1 ) );
  833.     new( rv1, makematrix( 1, 1 ) );
  834.     
  835.     If A^.r < A^.c Then aa := matequals( aa, Tran( A ) )
  836.     Else aa := matequals( aa, A );
  837.     
  838.     uu := matequals( uu, Fill( aa^.r, aa^.r, 0.0 ) );
  839.     vv := matequals( vv, Fill( aa^.c, aa^.c, 0.0 ) );
  840.     ss := matequals( ss, Fill( aa^.c, 1, 0.0 ) );
  841.     rv1 := matequals( rv1, Fill( aa^.c, 1, 0.0 ) );
  842.     
  843.     svdtemp( aa, ss, makeu, uu, makev, vv, ierr, rv1 );
  844.     
  845.     If A^.r < A^.c Then Begin 
  846.         If makeu Then U := matequals( u, vv );
  847.         If makev Then V := matequals( v, uu );
  848.     End 
  849.     Else Begin 
  850.         If makeu Then U := matequals( u, uu );
  851.         If makev Then V := matequals( v, vv );
  852.     End;
  853.     
  854.     S := matequals( s, ss );
  855.     dispose( aa, killvmatrix );
  856.     dispose( uu, killvmatrix );
  857.     dispose( ss, killvmatrix );
  858.     dispose( vv, killvmatrix );
  859.     dispose( rv1, killvmatrix );
  860.     
  861.     Dispatch^.Declevel;
  862.     svd := ierr;
  863. End;
  864.  
  865. { ---------------- end svd ------------------------------------  }
  866.  
  867. Function Ginv{(ROp : vmatrixptr): vmatrixptr;};
  868. Var
  869.   u, s, v, g : vmatrixptr;
  870.   i : integer;
  871.   t : double;
  872. Begin
  873.     ROp^.Garbage;
  874.     Dispatch^.Inclevel;
  875.     
  876.     new( u, makematrix( 1, 1 ) );
  877.     new( s, makematrix( 1, 1 ) );
  878.     new( v, makematrix( 1, 1 ) );
  879.     new( g, makematrix( 1, 1 ) );
  880.     
  881.     Svd( ROp, u, s, v, TRUE, TRUE );
  882.     For i := 1 To  s^.r Do Begin 
  883.         t := s^.m( i, 1 );
  884.         s^.mm( i, 1 )^ := 0;
  885.         If abs( t ) <> 0.0 Then s^.mm( i, 1 )^ := 1.0 / t;
  886.     End;
  887.     g := matequals( g, mult( mult( v, Diag( s ) ), Tran( u ) ) );
  888.     
  889.     dispose( u, killvmatrix );
  890.     dispose( s, killvmatrix );
  891.     dispose( v, killvmatrix );
  892.     Dispatch^.Push( g );
  893.     Ginv := Dispatch^.DecReturn;
  894. End;
  895.  
  896.  
  897. {  / -------------------------- fft ------------------------------  }
  898. {  /  de Boor (1980) SIAM J SCI. STAT. COMPUT. V1 no 1, pp 173-178  }
  899. {  /  and NewMat by R. G. Davies a C++ matrix Package  }
  900. {  / -------------------------------------------------------------  }
  901.  
  902. Procedure cossin( n, d: integer; Var c, s: double );
  903. Var
  904. { calculate cos(twopi*n/d) and sin(twopi*n/d)  }
  905. { minimise roundoff error  }
  906.         n4 : longint;
  907.         sector, sign : integer;
  908.         ratio, dn4, dd, divis  : double;
  909. Begin
  910.     n4 := n * 4;
  911.     dn4 := n4;
  912.     dd := d;
  913.     divis := dn4 / dd;
  914.     If divis < 0 Then divis := divis - 0.5 
  915.     Else divis := divis + 0.5;
  916.     sector := trunc( divis );
  917.     n4 := n4 - sector * d;
  918.     dn4 := n4;
  919.     If sector < 0 Then sector := 3 - (3 - sector) Mod 4
  920.     Else sector := sector Div 4;
  921.     ratio := 1.5707963267948966192 * dn4 / dd;
  922.     
  923.     Case sector Of
  924.         0 : Begin 
  925.             c := cos( ratio );   s := sin( ratio );  
  926.         End;
  927.         1 : Begin 
  928.             c := - sin( ratio );  s := cos( ratio );  
  929.         End;
  930.         2 : Begin 
  931.             c := - cos( ratio );  s := - sin( ratio ); 
  932.         End;
  933.         3 : Begin 
  934.             c := sin( ratio );   s := - cos( ratio ); 
  935.         End;
  936.         Else writeln( 'got here' );
  937.     End;
  938. End;
  939.  
  940. Procedure fftstep( Var AB, XY: vmatrixptr; after, now, before: integer );
  941. Var
  942.   gamma, delta, x, m, j, a, ia, x1, a1, x2, ib, a2, iin : integer;
  943.   r_arg, i_arg, r_value, i_value, temp : double;
  944. Begin
  945.     gamma := after * before;
  946.     delta := now * after;
  947.     
  948.     r_arg := 1.0; i_arg := 0.0;
  949.     
  950.     x := 1;
  951.     m := AB^.r - gamma;
  952.     
  953.     For j := 0 To now - 1 Do Begin 
  954.         a := 1;
  955.         x1 := x;
  956.         x := x + after;
  957.         For ia := 0 To after - 1 Do Begin 
  958.             {   /  generate sins & cosines explicitly rather than iteratively  }
  959.             {   /  for more accuracy; but slower  }
  960.             cossin( - (j * after + ia), delta, r_arg, i_arg );
  961.             
  962.             a1 := a;  a := a + 1;
  963.             x2 := x1; x1 := x1 + 1;
  964.             If now = 2 Then Begin 
  965.                 ib := before;
  966.                 While ib <> 0 Do Begin 
  967.                     ib := ib - 1;
  968.                     a2 := m + a1;
  969.                     a1 := a1 + after;
  970.                     r_value := AB^.m( a2, 1 );
  971.                     i_value := AB^.m( a2, 2 );
  972.                     XY^.mm( x2, 1 )^ := r_value * r_arg - i_value * i_arg + AB^.m( a2 - gamma, 1 );
  973.                     XY^.mm( x2, 2 )^ := r_value * i_arg + i_value * r_arg + AB^.m( a2 - gamma, 2 );
  974.                     x2 := x2 + delta;
  975.                 End;
  976.             End
  977.             Else Begin 
  978.                 ib := before;
  979.                 While ib <> 0 Do Begin 
  980.                     ib := ib - 1;
  981.                     a2 := m + a1;
  982.                     a1 := a1 + after;
  983.                     r_value := AB^.m( a2, 1 );
  984.                     i_value := AB^.m( a2, 2 );
  985.                     iin := now - 1;
  986.                     While iin <> 0 Do Begin 
  987.                         iin := iin - 1;
  988.                         a2 := a2 - gamma;
  989.                         temp := r_value;
  990.                         r_value := r_value * r_arg - i_value * i_arg + AB^.m( a2, 1 );
  991.                         i_value := temp * i_arg + i_value * r_arg + AB^.m( a2, 2 );
  992.                     End;
  993.                     XY^.mm( x2, 1 )^ := r_value;
  994.                     XY^.mm( x2, 2 )^ := i_value;
  995.                     x2 := x2 + delta;
  996.                 End;
  997.             End;
  998.         End;
  999.     End;
  1000. End;
  1001.  
  1002.  
  1003. Function Fft{(ROp: vmatrixptr; INZEE : boolean): vmatrixptr;};
  1004. {  /  if INZEE = TTRUE  then calculate fft  }
  1005. {  /  if INZEE = FFALSE then calculate inverse fft  }
  1006. Var
  1007.   n, nextmx, after, before, next, b1, now, i: integer;
  1008.   prime : Array[1..26] Of integer;
  1009.   dn : double;
  1010.   AB, XY : vmatrixptr;
  1011.   iinzee : boolean;
  1012. Label foundit;
  1013. Begin
  1014.     ROp^.Garbage;
  1015.     If (ROp^.c <> 1) And (ROp^.c <> 2 ) Then
  1016.         Dispatch^.Nrerror( 'Fft: Input must have 1 or 2 columns' );
  1017.     
  1018.     n := ROp^.r;                       {// length of arrays }
  1019.     dn := n;
  1020.     nextmx := 26;
  1021.     
  1022.     prime[1] := 2;   prime[2] := 3;   prime[3] := 5;   prime[4] := 7;
  1023.     prime[5] := 11;  prime[6] := 13;  prime[7] := 17;  prime[8] := 19;
  1024.     prime[9] := 23;  prime[10] := 29; prime[11] := 31; prime[12] := 37;
  1025.     prime[13] := 41; prime[14] := 43; prime[15] := 47; prime[16] := 53;
  1026.     prime[17] := 59; prime[18] := 61; prime[19] := 67; prime[20] := 71;
  1027.     prime[21] := 73; prime[22] := 79; prime[23] := 83; prime[24] := 89;
  1028.     prime[25] := 97; prime[26] := 101;
  1029.     
  1030.     after := 1;
  1031.     before := n;
  1032.     next := 1;
  1033.     iinzee := TRUE;
  1034.     
  1035.     Dispatch^.Inclevel;
  1036.     new( AB, makematrix( n, 2 ) );
  1037.     new( XY, makematrix( n, 2 ) );
  1038.     
  1039.     If ROp^.c = 1 Then Begin 
  1040.         For i := 1 To  n Do Begin 
  1041.             AB^.mm( i, 1 )^ := ROp^.m( i, 1 );
  1042.             AB^.mm( i, 2 )^ := 0.0;
  1043.         End;
  1044.     End
  1045.     Else AB := matequals( AB, ROp );
  1046.     XY := matequals( XY, Fill( n, 2, 0.0 ) );
  1047.     
  1048.     If Not INZEE Then Begin 
  1049.         {  take complex conjugate for ifft  }
  1050.         For i := 1 To n Do AB^.mm( i, 2 )^ := - AB^.m( i, 2 );
  1051.     End;
  1052.     
  1053.     Repeat
  1054.         While TRUE Do Begin 
  1055.             If next <= nextmx Then now := prime[next];
  1056.             b1 := before Div now;
  1057.             If b1 * now = before Then Goto foundit;
  1058.             next := next + 1;
  1059.             now := now + 2;
  1060.         End;
  1061.         foundit : before := b1;
  1062.         
  1063.         If iinzee = TRUE Then fftstep( AB, XY, after, now, before )
  1064.         Else fftstep( XY, AB, after, now, before );
  1065.         
  1066.         {iinzee :=(iinzee = TTRUE) ? FFALSE : TTRUE;}
  1067.         If iinzee = TRUE Then iinzee := FALSE 
  1068.         Else iinzee := TRUE;
  1069.         after := after * now;
  1070.     Until before = 1 ;
  1071.     
  1072.     If iinzee = TRUE Then XY := matequals( XY, AB );
  1073.     {  divide by n for ifft  }
  1074.     If Not INZEE Then
  1075.         For i := 1 To XY^.r Do Begin 
  1076.             XY^.mm( i, 1 )^ := XY^.m( i, 1 ) / dn;
  1077.             XY^.mm( i, 2 )^ := XY^.m( i, 2 ) / dn;
  1078.         End;
  1079.     
  1080.     Dispatch^.Push( XY );
  1081.     dispose( AB, killvmatrix );
  1082.     fft := Dispatch^.DecReturn;
  1083. End;
  1084.  
  1085.  
  1086.  
  1087.  
  1088. {  ///////////////////   reshaping functions  }
  1089.  
  1090. Function Vec{( ROp : vmatrixptr): vmatrixptr;};
  1091. Var
  1092.   lrc, lr, lc : longint;
  1093.   r,c,cnter, i,j : integer;
  1094.   temp : vmatrixptr;
  1095. Begin
  1096.     {// send columns of ROp to a vector}
  1097.     ROp^.Garbage;
  1098.     
  1099.     lr := r;
  1100.     lc := c;
  1101.     lrc := ROp^.r * ROp^.c;
  1102.     If (lrc >= 32768) Or (lrc < 1) Then
  1103.         Dispatch^.Nrerror( ' Vec: vec produces invalid indices' );
  1104.     r := ROp^.r * ROp^.c;
  1105.     c := 1;
  1106.     new( temp, makematrix( r, c ) );
  1107.     cnter := 1;
  1108.     For j := 1 To ROp^.c Do
  1109.         For i := 1 To ROp^.r Do Begin 
  1110.             temp^.mm( cnter, 1 )^ := ROp^.m( i, j );
  1111.             cnter := cnter + 1;
  1112.         End;
  1113.     
  1114.     Dispatch^.Push( temp );
  1115.     vec := Dispatch^.ReturnMat;
  1116. End;
  1117.  
  1118. Function Diag{(ROp : vmatrixptr) : vmatrixptr;};
  1119. Var
  1120.   temp : vmatrixptr;
  1121.   i, rr, cc : integer;
  1122. Begin
  1123.     { make a diagonal matrix from a vector or the principal diagonal of}
  1124.     { a square matrix, off diagonal elements are zero  }
  1125.     ROp^.Garbage;
  1126.     If (ROp^.r <> ROp^.c) And  (ROp^.c <> 1) Then
  1127.         Dispatch^.Nrerror( 'Diag: ROp is not square or is not a column vector' );
  1128.     Dispatch^.Inclevel;
  1129.     
  1130.     new( temp, makematrix( 1, 1 ) );
  1131.     temp := matequals( temp, Fill( ROp^.r, ROp^.r, 0.0 ) );
  1132.     
  1133.     rr := ROp^.r;
  1134.     cc := ROp^.c;
  1135.     For i := 1 To rr Do
  1136.         If rr = cc Then temp^.mm( i, i )^ := ROp^.m( i, i )
  1137.         Else temp^.mm( i, i )^ := ROp^.m( i, 1 );
  1138.     
  1139.     Dispatch^.Push( temp );
  1140.     Diag := Dispatch^.DecReturn;
  1141. End;
  1142.  
  1143.  
  1144. Function Shape{(ROp: vmatrixptr; nrows: integer) : vmatrixptr;};
  1145. Var
  1146.  nr, lr, lc, lrc : longint;
  1147.  cnter, i,j : integer;
  1148.  temp1, temp: vmatrixptr;
  1149. Begin
  1150.     {reshape a matrix into n rows, nrows must divide r*c without a}
  1151.     {remainder  }
  1152.     ROp^.Garbage;
  1153.     
  1154.     nr := nrows;
  1155.     lr := ROp^.r;
  1156.     lc := ROp^.c;
  1157.     
  1158.     If (lr * lc Mod nr) <> 0 Then
  1159.         Dispatch^.Nrerror( ' Shape: nrows divides r*c with a remainder' );
  1160.     
  1161.     lrc := (lr * lc) Div nr;
  1162.     If (nr >= 32768) Or (nr < 1) Or (lrc >= 32768) Or (lrc < 1) Then
  1163.         Dispatch^.Nrerror( ' Shape: reshape produces invalid indices' );
  1164.     
  1165.     Dispatch^.Inclevel;
  1166.     new( temp1, makematrix( 1, 1 ) );
  1167.     temp1 := matequals( temp1, Vec( ROp ) );
  1168.     
  1169.     new( temp, makematrix( nrows, lrc ) );
  1170.     cnter := 1;
  1171.     For i := 1 To temp^.r Do
  1172.         For j := 1 To temp^.c Do Begin 
  1173.             temp^.mm( i, j )^ := temp1^.m( cnter, 1 );
  1174.             cnter := cnter + 1;
  1175.         End;
  1176.     
  1177.     Dispatch^.Push( temp );
  1178.     dispose( temp1, killvmatrix );
  1179.     Shape := Dispatch^.DecReturn;
  1180. End;
  1181.  
  1182.  
  1183. Function Sum{( ROp : vmatrixptr; marg : margins ) : vmatrixptr;};
  1184. Var
  1185.   i,j,r,c : integer;
  1186.   ssum : double;
  1187.   temp : vmatrixptr;
  1188. Begin
  1189.     { sum(ROp, ALL) returns 1x1}
  1190.     { sum(ROp, ROWS) returns 1xc  }
  1191.     { sum(ROp, COLUMNS) returns rx1  }
  1192.     ROP^.Garbage;
  1193.     
  1194.     Case marg Of
  1195.         ALL     : Begin 
  1196.             r := 1; c := 1; 
  1197.         End;
  1198.         ROWS    : Begin 
  1199.             r := 1; c := ROP^.c; 
  1200.         End;
  1201.         COLUMNS : Begin 
  1202.             r := ROP^.r; c := 1; 
  1203.         End;
  1204.         Else
  1205.             Dispatch^.Nrerror( ' Sum: invalid margin specification' );
  1206.     End;
  1207.     new( temp, makematrix( r, c ) );
  1208.     
  1209.     Case marg Of
  1210.         ALL : Begin 
  1211.             ssum := 0.0;
  1212.             For i := 1 To ROp^.r Do
  1213.                 For j := 1 To ROp^.c Do ssum := ssum + ROp^.m( i, j );
  1214.             temp^.mm( 1, 1 )^ := ssum;
  1215.         End;
  1216.         ROWS : Begin  
  1217.             For j := 1 To ROp^.c Do Begin 
  1218.                 ssum := 0.0;
  1219.                 For i := 1 To ROp^.r Do ssum := ssum + ROp^.m( i, j );
  1220.                 temp^.mm( 1, j )^ := ssum;
  1221.             End;
  1222.         End;
  1223.         COLUMNS : Begin 
  1224.             For i := 1 To ROp^.r Do Begin 
  1225.                 ssum := 0.0;
  1226.                 For j := 1 To ROp^.c Do ssum := ssum + ROp^.m( i, j );
  1227.                 temp^.mm( i, 1 )^ := ssum;
  1228.             End;
  1229.         End;
  1230.     End;
  1231.     
  1232.     Dispatch^.Push( temp );
  1233.     Sum := Dispatch^.ReturnMat;
  1234. End;
  1235.  
  1236. Function Sumsq{(ROp: vmatrixptr; marg: margins): vmatrixptr;};
  1237. Var
  1238.   i,j,r,c : integer;
  1239.   sum : double;
  1240.   temp : vmatrixptr;
  1241. Begin
  1242.     { sumsq(ROp, ALL) returns 1x1}
  1243.     { sumsq(ROp, ROWS) returns 1xc  }
  1244.     { sumsq(ROp, COLUMNS) returns cx1  }
  1245.     ROp^.Garbage;
  1246.     
  1247.     Case marg Of
  1248.         ALL : Begin 
  1249.             r := 1; c := 1; 
  1250.         End;
  1251.         ROWS : Begin 
  1252.             r := 1; c := ROp^.c; 
  1253.         End;
  1254.         COLUMNS : Begin 
  1255.             r := ROp^.r; c := 1; 
  1256.         End;
  1257.         Else
  1258.             Dispatch^.Nrerror( ' Sumsq: invalid margin specification' );
  1259.     End;
  1260.     new( temp, makematrix( r, c ) );
  1261.     
  1262.     Case marg Of
  1263.         ALL : Begin 
  1264.             sum := 0.0;
  1265.             For i := 1 To ROp^.r Do
  1266.                 For j := 1 To ROp^.c Do sum := sum + sqr( ROp^.m( i, j ) );
  1267.             temp^.mm( 1, 1 )^ := sum;
  1268.         End;
  1269.         ROWS : Begin 
  1270.             For j := 1 To  ROp^.c Do Begin 
  1271.                 sum := 0.0;
  1272.                 For i := 1 To ROp^.r Do
  1273.                     sum := sum + sqr( ROp^.m( i, j ) );
  1274.                 temp^.mm( 1, j )^ := sum;
  1275.             End;
  1276.         End;
  1277.         COLUMNS : Begin 
  1278.             For i := 1 To ROp^.r Do Begin    
  1279.                 sum := 0.0;
  1280.                 For j := 1 To ROp^.c Do sum := sum + sqr( ROp^.m( i, j ) );
  1281.                 temp^.mm( i, 1 )^ := sum;
  1282.             End;
  1283.         End;
  1284.     End;
  1285.     
  1286.     Dispatch^.Push( temp );
  1287.     sumsq := Dispatch^.ReturnMat;
  1288. End;
  1289.  
  1290. Function Cusum{( ROp: vmatrixptr): vmatrixptr;};
  1291. Var
  1292.   i,j : integer;
  1293.   sum : double;
  1294.   temp : vmatrixptr;
  1295. Begin
  1296.     {// cumulative sum along rows}
  1297.     ROp^.Garbage;
  1298.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  1299.     
  1300.     sum := 0.0;
  1301.     For i := 1 To ROp^.r Do Begin 
  1302.         For j := 1 To ROp^.c Do Begin 
  1303.             sum := sum + ROp^.m( i, j );
  1304.             temp^.mm( i, j )^ := sum;
  1305.         End;
  1306.     End;
  1307.     Dispatch^.Push( temp );
  1308.     Cusum := Dispatch^.ReturnMat;
  1309. End;
  1310.  
  1311. Function Mmax{( ROp: vmatrixptr; marg : margins): vmatrixptr;};
  1312. Var
  1313.   i,j,r,c,maxr,maxc : integer;
  1314.   mmmax : double;
  1315.   temp : vmatrixptr;
  1316. Begin
  1317.     {// Mmax(ROp, ALL) returns 3x1}
  1318.     {/  Mmax(ROp, ROWS) returns 3xc  }
  1319.     {/  Mmax(ROp, COLUMNS) returns rx3  }
  1320.     ROp^.Garbage;
  1321.     
  1322.     Case marg Of
  1323.         ALL : Begin     
  1324.             r := 3; c := 1; 
  1325.         End;
  1326.         ROWS : Begin    
  1327.             r := 3; c := ROp^.c; 
  1328.         End;
  1329.         COLUMNS : Begin 
  1330.             r := ROp^.r; c := 3; 
  1331.         End;
  1332.         Else
  1333.             Dispatch^.Nrerror( ' Mmax: invalid margin specification' );
  1334.     End;
  1335.     new( temp, makematrix( r, c ) );
  1336.     
  1337.     Case marg Of
  1338.         ALL : Begin 
  1339.             mmmax := ROp^.m( 1, 1 );
  1340.             maxr := 1;
  1341.             maxc := 1;
  1342.             For i := 1 To  ROp^.r Do
  1343.                 For j := 1 To  ROp^.c Do
  1344.                     If mmmax < ROp^.m( i, j ) Then Begin
  1345.                         mmmax := ROp^.m( i, j );
  1346.                         maxr := i;
  1347.                         maxc := j;
  1348.                     End;
  1349.             temp^.mm( 1, 1 )^ := maxr;
  1350.             temp^.mm( 2, 1 )^ := maxc;
  1351.             temp^.mm( 3, 1 )^ := mmmax;
  1352.         End;
  1353.         ROWS : Begin 
  1354.             For j := 1 To c Do Begin 
  1355.                 mmmax := ROp^.m( 1, j );
  1356.                 maxr := 1;
  1357.                 maxc := j;
  1358.                 For i := 1 To ROp^.r Do
  1359.                     If mmmax < ROp^.m( i, j ) Then Begin
  1360.                         maxr := i;
  1361.                         mmmax := ROp^.m( i, j );
  1362.                     End;
  1363.                 temp^.mm( 1, j )^ := maxr;
  1364.                 temp^.mm( 2, j )^ := maxc;
  1365.                 temp^.mm( 3, j )^ := mmmax;
  1366.             End;
  1367.         End;
  1368.         COLUMNS : Begin 
  1369.             For i := 1 To r Do Begin 
  1370.                 mmmax := ROp^.m( i, 1 );
  1371.                 maxr := i;
  1372.                 maxc := 1;
  1373.                 For j := 1 To ROp^.c Do
  1374.                     If mmmax < ROp^.m( i, j ) Then Begin 
  1375.                         maxc := j;
  1376.                         mmmax := ROp^.m( i, j );
  1377.                     End;
  1378.                 temp^.mm( i, 1 )^ := maxr;
  1379.                 temp^.mm( i, 2 )^ := maxc;
  1380.                 temp^.mm( i, 3 )^ := mmmax;
  1381.             End;
  1382.         End;
  1383.     End;
  1384.     Dispatch^.Push( temp );
  1385.     Mmax := Dispatch^.ReturnMat;
  1386. End;
  1387.  
  1388. Function Mmin{( ROp: vmatrixptr; marg : margins): vmatrixptr;};
  1389. Var
  1390.   i,j,r,c,minr,minc : integer;
  1391.   mmmin : double;
  1392.   temp : vmatrixptr;
  1393. Begin
  1394.     {/  Mmin(ROp, ALL) returns 3x1}
  1395.     {/  Mmin(ROp, ROWS) returns 3xc  }
  1396.     {/  Mmin(ROp, COLUMNS) returns rx3  }
  1397.     ROp^.Garbage;
  1398.     
  1399.     Case marg Of
  1400.         ALL : Begin     
  1401.             r := 3; c := 1; 
  1402.         End;
  1403.         ROWS : Begin    
  1404.             r := 3; c := ROp^.c; 
  1405.         End;
  1406.         COLUMNS : Begin 
  1407.             r := ROp^.r; c := 3; 
  1408.         End;
  1409.         Else
  1410.             Dispatch^.Nrerror( ' Mmin: invalid margin specification' );
  1411.     End;
  1412.     new( temp, makematrix( r, c ) );
  1413.     
  1414.     Case marg Of
  1415.         ALL : Begin 
  1416.             mmmin := ROp^.m( 1, 1 );
  1417.             minr := 1;
  1418.             minc := 1;
  1419.             For i := 1 To  ROp^.r Do
  1420.                 For j := 1 To  ROp^.c Do
  1421.                     If mmmin > ROp^.m( i, j ) Then Begin
  1422.                         mmmin := ROp^.m( i, j );
  1423.                         minr := i;
  1424.                         minc := j;
  1425.                     End;
  1426.             temp^.mm( 1, 1 )^ := minr;
  1427.             temp^.mm( 2, 1 )^ := minc;
  1428.             temp^.mm( 3, 1 )^ := mmmin;
  1429.         End;
  1430.         ROWS : Begin 
  1431.             For j := 1 To c Do Begin 
  1432.                 mmmin := ROp^.m( 1, j );
  1433.                 minr := 1;
  1434.                 minc := j;
  1435.                 For i := 1 To ROp^.r Do
  1436.                     If mmmin > ROp^.m( i, j ) Then Begin
  1437.                         minr := i;
  1438.                         mmmin := ROp^.m( i, j );
  1439.                     End;
  1440.                 temp^.mm( 1, j )^ := minr;
  1441.                 temp^.mm( 2, j )^ := minc;
  1442.                 temp^.mm( 3, j )^ := mmmin;
  1443.             End;
  1444.         End;
  1445.         COLUMNS : Begin 
  1446.             For i := 1 To r Do Begin 
  1447.                 mmmin := ROp^.m( i, 1 );
  1448.                 minr := i;
  1449.                 minc := 1;
  1450.                 For j := 1 To ROp^.c Do
  1451.                     If mmmin > ROp^.m( i, j ) Then Begin 
  1452.                         minc := j;
  1453.                         mmmin := ROp^.m( i, j );
  1454.                     End;
  1455.                 temp^.mm( i, 1 )^ := minr;
  1456.                 temp^.mm( i, 2 )^ := minc;
  1457.                 temp^.mm( i, 3 )^ := mmmin;
  1458.             End;
  1459.         End;
  1460.     End;
  1461.     Dispatch^.Push( temp );
  1462.     Mmin := Dispatch^.ReturnMat;
  1463. End;
  1464.  
  1465. {  ///////////////////////////////  elemenatary row and column operations  }
  1466. {  /////////////////  note these functions operate directly on ROp  }
  1467.  
  1468. Procedure Crow{(var ROp: vmatrixptr; row : integer; c : double)};
  1469. Var
  1470.   i : integer;
  1471. Begin
  1472.     {// multiply a row by a constant}
  1473.     ROp^.Garbage;
  1474.     
  1475.     If (row < 1) Or (row > ROp^.r) Then
  1476.         Dispatch^.Nrerror( ' Crow: row index out of range' );
  1477.     
  1478.     For i := 1 To ROp^.c Do ROp^.mm( row, i )^ := c * ROp^.m( row, i );
  1479. End;
  1480.  
  1481. Procedure Srow{(var ROp : vmatrixptr; row1, row2 : integer);};
  1482. Var
  1483.   i,j : integer;
  1484.   tmp : double;
  1485. Begin
  1486.     {// swap rows }
  1487.     ROp^.Garbage;
  1488.     If (row1 < 1) Or (row1 > ROp^.r) Or (row2 < 1) Or (row2 > ROp^.r) Then Dispatch^.Nrerror( ' Srow: row index out of range' );
  1489.     
  1490.     If row1 = row2 Then exit;
  1491.     For i := 1 To ROp^.c Do Begin 
  1492.         tmp := ROp^.m( row1, i );
  1493.         ROp^.mm( row1, i )^ := ROp^.m( row2, i );
  1494.         ROp^.mm( row2, i )^ := tmp;
  1495.     End;
  1496. End;
  1497.  
  1498. Procedure Lrow{(var ROp : vmatrixptr; row1, row2 : integer; c : double);};
  1499. Var
  1500.   i : integer;
  1501. Begin
  1502.     {// add c*r1 to r2}
  1503.     ROp^.Garbage;
  1504.     If (row1 < 1) Or (row1 > ROp^.r) Or (row2 < 1) Or (row2 > ROp^.r) Then Dispatch^.Nrerror( ' Lrow: row index out of range' );
  1505.     
  1506.     For i := 1 To ROp^.c Do
  1507.         ROp^.mm( row2, i )^ := ROp^.m( row2, i ) + c * ROp^.m( row1, i );
  1508. End;
  1509.  
  1510. Procedure Ccol{(var ROp : vmatrixptr; col : integer; c : double);};
  1511. Var
  1512.   i : integer;
  1513. Begin
  1514.     {// multiply a col by a constant}
  1515.     ROp^.Garbage;
  1516.     
  1517.     If (col < 1) Or (col > ROp^.c) Then
  1518.         Dispatch^.Nrerror( ' Ccol: col index out of range' );
  1519.     
  1520.     For i := 1 To  ROp^.r Do ROp^.mm( i, col )^ := c * ROp^.m( i, col );
  1521. End;
  1522.  
  1523. Procedure Scol{(var ROp : vmatrixptr; col1, col2 : integer);};
  1524. Var
  1525.   i : integer;
  1526.   tmp : double;
  1527. Begin
  1528.     {// swap cols}
  1529.     ROp^.Garbage;
  1530.     If (col1 < 1) Or (col1 > ROp^.c) Or (col2 < 1) Or (col2 > ROp^.c) Then Dispatch^.Nrerror( ' Scol: col index out of range' );
  1531.     
  1532.     If col1 = col2 Then exit;
  1533.     For i := 1 To ROp^.r Do Begin 
  1534.         tmp := ROp^.m( i, col1 );
  1535.         ROp^.mm( i, col1 )^ := ROp^.m( i, col2 );
  1536.         ROp^.mm( i, col2 )^ := tmp;
  1537.     End;
  1538. End;
  1539.  
  1540. Procedure Lcol{(var ROp : vmatrixptr; col1, col2 : integer; c : double);};
  1541. Var
  1542.   i: integer;
  1543. Begin
  1544.     {// add c*c1 to c2}
  1545.     ROp^.Garbage;
  1546.     If (col1 < 1) Or (col1 > ROp^.c) Or (col2 < 1) Or (col2 > ROp^.c) Then Dispatch^.Nrerror( ' Lcol: col index out of range' );
  1547.     
  1548.     For i := 1 To ROp^.r Do
  1549.         ROp^.mm( i, col2 )^ := ROp^.m( i, col2 ) + c * ROp^.m( i, col1 );
  1550. End;
  1551.  
  1552. Begin
  1553. End.
  1554.  
  1555.